Sys.setlocale('LC_ALL', 'C')
[1] "LC_CTYPE=C;LC_NUMERIC=C;LC_TIME=C;LC_COLLATE=C;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=vi_VN;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=vi_VN;LC_IDENTIFICATION=C"
library(fpp2)
library(forecast)
library(readxl)
library(DIMORA)
setwd('/home/hieu/BDMA/TSA')
Competition between cassettes and CDs
music<- read_excel('music.xlsx')
colors <- grey(seq(0, 1, length = 5))
year<- music$year[1:36]
cass<- music$cassette[1:36]
cd<- music$cd[10:36]
#
##Exploratory plots
#instantaneous
plot(year[10:36],cass[10:36],xlab='Year', ylab='Annual sales', type= 'b', pch=16, lty=3, cex=0.7, ylim=c(0,1000), col=colors[1])
points(year[10:36],cd, type= 'b', pch=2, lty=3, cex=0.7, col=colors[3])
legend('topleft', legend=c('Cassette', 'CD'), pch=c(16,2), col=c(colors[1], colors[3]))
#
#cumulative
plot(year[10:36],sum(cass[1:9])+cumsum(cass[10:36]),xlab='Year', ylab='Cumulative sales', type= 'b', pch=16, lty=3, cex=0.7, ylim=c(0,15000), col=colors[1])
points(year[10:36],cumsum(cd), type= 'b', pch=2, lty=3, cex=0.7, col=colors[3])
legend('topleft', legend=c('Cassette', 'CD'), pch=c(16,2), col=c(colors[1], colors[3]))
# UCRCD model with delta and gamma (par = 'double')
ucrcd<- UCRCD(cass,cd, display=T)
summary(ucrcd)
Call: ( UCRCD Model )
UCRCD(series1 = cass, series2 = cd, display = T)
Residuals Series 1:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-61.911 -5.930 2.142 1.660 14.593 55.344
Residuals Series 2:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-84.0608 -13.5927 0.8525 0.1994 18.4845 57.6658
Coefficients:
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error Series 1: 28.17771 on 26 degrees of freedom
Residual standard error Series 2: 43.43298 on 17 degrees of freedom
Multiple R-squared: 0.989805 Residual sum of squares: 52712.77
# Plot observed vs fitted
plot(year[10:36],cass[10:36],xlab='Year', ylab='Annual sales', type= 'b', pch=16, lty=3, cex=0.7, ylim=c(0,1000), col=colors[1])
points(year[10:36],cd, type= 'b', pch=2, lty=3, cex=0.7, col=colors[3])
lines(year[10:36], (ucrcd$fitted.i[[1]])[10:36], lwd=1, col=colors[1])
lines(year[10:36], (ucrcd$fitted.i[[2]]), lwd=1, col=colors[3])
legend('topleft', legend=c('Cassette', 'CD'), lty=1, lwd=2, col=colors[c(1,3)], pch=c(16,2), cex=0.7)
plot(year[10:36],sum(cass[1:9])+cumsum(cass[10:36]),xlab='Year', ylab='Cumulative sales', type= 'b', pch=16, lty=3, cex=0.7, ylim=c(0,15000), col=colors[1])
points(year[10:36],cumsum(cd), type= 'b', pch=2, lty=3, cex=0.7, col=colors[3])
lines(year[10:36], cumsum(ucrcd$fitted.i[[1]])[10:36], lwd=1, col=colors[1])
lines(year[10:36], cumsum(ucrcd$fitted.i[[2]]), lwd=1, col=colors[3])
legend('topleft', legend=c('Cassette', 'CD'), lty=1, lwd=2, col=colors[c(1,3)], pch=c(16,2), cex=0.7)
Data on US personal consumption and income
uschange
Consumption Income Production Savings Unemployment
1970 Q1 0.61598622 0.97226104 -2.45270031 4.81031150 0.9
1970 Q2 0.46037569 1.16908472 -0.55152509 7.28799234 0.5
1970 Q3 0.87679142 1.55327055 -0.35870786 7.28901306 0.5
1970 Q4 -0.27424514 -0.25527238 -2.18545486 0.98522964 0.7
1971 Q1 1.89737076 1.98715363 1.90973412 3.65777061 -0.1
1971 Q2 0.91199291 1.44733417 0.90153584 6.05134180 -0.1
1971 Q3 0.79453885 0.53181193 0.30801942 -0.44583221 0.1
1971 Q4 1.64858747 1.16012514 2.29130441 -1.53087186 0.0
1972 Q1 1.31372218 0.45701150 4.14957387 -4.35859438 -0.2
1972 Q2 1.89147495 1.01662441 1.89062398 -5.05452579 -0.1
1972 Q3 1.53071400 1.90410126 1.27335290 5.80995904 -0.2
1972 Q4 2.31829471 3.89025866 3.43689207 16.04471706 -0.3
1973 Q1 1.81073916 0.70825266 2.79907636 -5.34886849 -0.3
1973 Q2 -0.04173996 0.79430954 0.81768862 8.42603436 0.0
1973 Q3 0.35423556 0.43381827 0.86899693 2.75879565 -0.1
1973 Q4 -0.29163216 1.09380979 1.47296187 11.14642986 0.1
1974 Q1 -0.87702794 -1.66168482 -0.88248358 -2.53351449 0.2
1974 Q2 0.35113555 -0.93835321 0.07427919 -6.59264464 0.3
1974 Q3 0.40959770 0.09448779 -0.41314971 0.51717884 0.5
1974 Q4 -1.47580863 -0.12259599 -4.06411893 11.34339540 1.3
1975 Q1 0.83225762 -0.16369546 -6.85103912 -5.47619069 1.4
1975 Q2 1.65583461 4.53650956 -1.33129558 24.30960536 0.2
1975 Q3 1.41942029 -1.46376532 2.42435972 -17.65616104 -0.4
1975 Q4 1.05437932 0.76166351 2.16904208 0.64809041 -0.2
1976 Q1 1.97998024 1.16825761 3.02720471 -2.95006644 -0.6
1976 Q2 0.91391607 0.51729906 1.27881101 -1.47455755 0.0
1976 Q3 1.05532326 0.73370026 1.30386487 -0.06754475 0.0
1976 Q4 1.29889825 0.59458339 1.77537765 -3.57672239 0.2
1977 Q1 1.13637586 -0.03108003 2.05516067 -9.16055658 -0.4
1977 Q2 0.54994073 1.23808955 3.05838507 9.09050404 -0.2
1977 Q3 0.94985262 1.51880293 1.10308888 7.94495719 -0.4
1977 Q4 1.49599724 1.91456240 0.63346850 6.69627648 -0.4
1978 Q1 0.57549599 0.70266687 -0.29339056 2.92296383 -0.1
1978 Q2 2.11120960 0.98314132 3.94815264 -6.81114259 -0.4
1978 Q3 0.41796279 0.71992620 0.87114701 4.79207162 0.1
1978 Q4 0.79792710 0.78553605 1.78447991 2.37118400 0.0
1979 Q1 0.50584598 1.05755946 0.42594327 7.77418337 -0.2
1979 Q2 -0.05775339 -0.86765105 -0.20491944 -5.28634896 -0.1
1979 Q3 0.97730010 0.47100340 -0.29723637 -1.84549644 0.2
1979 Q4 0.26826982 0.44037974 0.33560928 4.04959810 0.1
1980 Q1 -0.15391875 0.33827686 0.41056141 5.86168864 0.3
1980 Q2 -2.27411019 -1.46388507 -4.30076832 8.24322919 1.3
1980 Q3 1.07188123 1.21301507 -1.64181977 5.70775044 -0.1
1980 Q4 1.31644941 1.94243865 3.78045520 9.15098787 -0.3
1981 Q1 0.52472770 -0.26813406 0.24627687 -5.68139002 0.2
1981 Q2 -0.01728203 -0.02363025 0.30977573 0.88183993 0.1
1981 Q3 0.40165150 2.02680183 0.91707444 15.99035721 0.1
1981 Q4 -0.75287620 0.19560628 -2.25457797 7.80550650 0.9
1982 Q1 0.65938376 0.11969888 -2.07131293 -3.34243955 0.5
1982 Q2 0.36854173 0.57548997 -1.24766384 2.19400166 0.6
1982 Q3 0.76954464 0.53484410 -1.40050430 0.03499563 0.5
1982 Q4 1.80876006 0.44938311 -1.90375664 -9.57651468 0.7
1983 Q1 0.96802954 0.85588425 1.14655720 0.34595460 -0.5
1983 Q2 1.95946831 0.70632719 2.17942248 -10.17004699 -0.2
1983 Q3 1.73949442 1.49810999 3.36771897 0.21217916 -0.9
1983 Q4 1.56389332 2.13138911 2.58168445 8.21600068 -0.9
1984 Q1 0.84526442 2.02348788 2.89709545 13.86918150 -0.5
1984 Q2 1.41504495 1.64921136 1.53821324 4.38900229 -0.6
1984 Q3 0.76546608 1.36163845 0.72128740 6.51686089 0.1
1984 Q4 1.31380062 0.81927319 0.04115557 -2.87544931 0.0
1985 Q1 1.68655320 -0.23895759 0.32353159 -18.71008389 -0.1
1985 Q2 0.93436990 1.90677905 0.07020996 11.82871950 0.2
1985 Q3 1.90256675 -0.33536283 -0.14046924 -23.57393474 -0.3
1985 Q4 0.25656565 1.14181151 0.57978813 11.36628338 -0.1
1986 Q1 0.84304279 1.23951110 0.58132135 5.86126836 0.2
1986 Q2 1.11177390 1.31938549 -0.57641778 3.27551734 0.0
1986 Q3 1.79499406 0.70477150 0.37249329 -10.09044542 -0.2
1986 Q4 0.63768446 0.17977925 1.13734778 -4.82920131 -0.4
1987 Q1 0.01569397 0.81973366 1.30758228 12.46424452 0.0
1987 Q2 1.37731686 -0.97505791 1.75000563 -29.52866718 -0.4
1987 Q3 1.15225712 1.80185055 1.84366200 12.32810406 -0.3
1987 Q4 0.21016439 1.32743427 2.40645058 16.63076101 -0.2
1988 Q1 1.76316026 1.44861875 0.92013121 -0.96896505 0.0
1988 Q2 0.73053714 1.02084894 0.87316353 5.67776867 -0.3
1988 Q3 0.85083233 0.95820336 0.38103668 3.64649867 0.0
1988 Q4 1.13789838 0.96207024 0.70292025 -0.19730358 -0.1
1989 Q1 0.46064152 1.22693023 0.43372685 10.01461545 -0.3
1989 Q2 0.46937808 -0.29489091 -0.36675732 -8.15576525 0.3
1989 Q3 0.98950145 0.67822897 -0.62142121 -2.48622554 0.0
1989 Q4 0.43942767 0.80025832 0.42443392 5.44681102 0.1
1990 Q1 0.85543417 0.83939484 0.68265169 2.87544931 -0.2
1990 Q2 0.31230451 0.59572848 0.77446547 5.10951644 0.0
1990 Q3 0.40261313 0.03740765 0.41944800 -3.17767248 0.7
1990 Q4 -0.75910716 -0.79479735 -1.57345296 -0.17953326 0.4
1991 Q1 -0.34535008 0.21183290 -1.91422028 6.49315257 0.5
1991 Q2 0.83564224 0.69043356 0.59131506 -0.30920615 0.1
1991 Q3 0.48439843 0.36205181 1.36255645 -0.14086493 0.0
1991 Q4 -0.02626579 0.85100324 0.21710308 11.34193010 0.4
1992 Q1 1.85996999 2.12421067 -0.13365365 7.23265150 0.1
1992 Q2 0.68354371 1.04095059 1.76874773 5.46708666 0.4
1992 Q3 1.07661214 0.43562041 0.76167388 -5.93646090 -0.2
1992 Q4 1.18372396 0.34210852 1.05024577 -5.88618856 -0.2
1993 Q1 0.37817936 0.55877186 0.87901471 2.63464703 -0.4
1993 Q2 0.89392729 0.17627103 0.21755108 -6.91664675 0.0
1993 Q3 1.09813766 0.05868803 0.40135891 -11.99337844 -0.3
1993 Q4 0.88122025 0.65496353 1.49618275 -1.83708870 -0.2
1994 Q1 1.14064791 0.69846579 1.22213656 -5.18600629 0.0
1994 Q2 0.77176225 1.05367166 1.78250275 5.15609751 -0.4
1994 Q3 0.77214364 0.59247377 1.26718100 -2.42215898 -0.2
1994 Q4 1.07014805 1.38110661 2.04370404 6.32351898 -0.4
1995 Q1 0.26420505 0.94873528 1.02552601 10.11514398 -0.1
1995 Q2 0.89311141 0.22780635 0.33785685 -10.60541172 0.2
1995 Q3 0.91264702 0.88957006 0.90043887 -0.11570727 0.0
1995 Q4 0.70025425 0.57591998 0.87467273 -2.90726686 0.0
1996 Q1 0.92360967 0.95255663 0.69285195 2.55933958 -0.1
1996 Q2 1.07997887 0.95161791 2.11134752 -0.75802112 -0.2
1996 Q3 0.60055799 0.79369738 1.24418680 3.33843952 -0.1
1996 Q4 0.78298122 0.52035746 1.35396890 -3.33843952 0.2
1997 Q1 1.04949253 0.99858552 1.86714700 0.61269338 -0.2
1997 Q2 0.45219855 0.85103564 1.48763922 6.17532322 -0.2
1997 Q3 1.69654264 1.18352222 2.28632066 -7.22796452 -0.1
1997 Q4 1.18062797 1.42325742 2.48091341 5.43456565 -0.2
1998 Q1 1.02693626 2.10753052 1.10343775 19.35335228 0.0
1998 Q2 1.75069399 1.38767133 0.65122238 -4.81709478 -0.2
1998 Q3 1.30596977 1.01464427 0.72551955 -3.12983982 0.1
1998 Q4 1.45888615 0.80893032 1.44421674 -9.14923404 -0.2
1999 Q1 0.94821191 0.89173174 1.10341663 1.88735718 -0.2
1999 Q2 1.46971415 0.24722185 0.98574261 -23.49652903 0.1
1999 Q3 1.12921436 0.66729226 0.90279881 -9.86264835 -0.1
1999 Q4 1.45748895 1.46092242 1.75533234 2.35825225 -0.2
2000 Q1 1.51106759 1.95061335 0.99682019 12.28684080 0.0
2000 Q2 0.95508878 1.03174349 1.23293805 1.28001748 0.0
2000 Q3 0.96797647 1.16178668 -0.10225268 2.57390229 -0.1
2000 Q4 0.88629738 0.33725343 -0.20388383 -13.16296208 0.0
2001 Q1 0.42159086 0.84865826 -1.35143911 13.22491995 0.4
2001 Q2 0.25689982 -0.08818148 -1.25954437 -6.89043916 0.2
2001 Q3 0.36381084 2.33678920 -1.44101744 41.66826457 0.5
2001 Q4 1.51630321 -1.24443353 -1.06013675 -56.75209674 0.7
2002 Q1 0.29958257 2.40331419 0.70916406 50.75796205 0.0
2002 Q2 0.50899032 0.50559877 1.54280957 0.87861837 0.1
2002 Q3 0.69667241 -0.12828194 0.59478143 -14.70397426 -0.1
2002 Q4 0.53634306 0.47941927 -0.05776556 1.58733492 0.3
2003 Q1 0.43826169 0.27834026 0.53922789 0.49744834 -0.1
2003 Q2 1.10719086 1.43729445 -0.69876172 7.00891625 0.4
2003 Q3 1.46377882 1.62544947 0.60727351 6.18413150 -0.2
2003 Q4 0.77334046 0.40353864 1.00599126 -6.89274778 -0.4
2004 Q1 0.96768535 0.72653162 0.65792806 -2.96152040 0.1
2004 Q2 0.64760607 0.98056746 0.57461780 8.30885627 -0.2
2004 Q3 0.95117167 0.52450113 0.56330030 -8.99318286 -0.2
2004 Q4 1.02041702 1.24238706 1.38522763 6.23585017 0.0
2005 Q1 0.76172556 -0.96827007 1.39435718 -42.28191228 -0.2
2005 Q2 1.08136588 0.78835467 0.50586367 -18.27592893 -0.2
2005 Q3 0.77186494 0.51136949 -0.50305848 -7.87665229 0.0
2005 Q4 0.37591485 0.82191843 0.93365010 20.37236078 -0.1
2006 Q1 1.11522822 2.25904474 0.95057853 37.40653542 -0.2
2006 Q2 0.53100554 0.14987813 0.59636010 -12.34810568 -0.1
2006 Q3 0.58208747 0.28490722 0.33552773 -10.55276140 -0.1
2006 Q4 1.01434389 1.30059162 0.25603401 6.03100080 -0.1
2007 Q1 0.52486184 0.65373993 0.91794957 6.60516929 0.0
2007 Q2 0.33874119 0.19260870 1.19594247 -7.23648452 0.2
2007 Q3 0.44391875 0.26238732 0.22356909 -9.00674555 0.1
2007 Q4 0.12505584 0.08392938 0.16424632 2.32887238 0.3
2008 Q1 -0.20652548 0.71926565 -0.42872571 29.83728599 0.1
2008 Q2 0.16783443 2.08693775 -1.41297022 46.43989041 0.5
2008 Q3 -0.72499446 -2.32611860 -3.26349945 -32.53252494 0.5
2008 Q4 -1.21068558 0.64019534 -4.35417741 36.31240490 1.2
2009 Q1 -0.34354370 -0.18888849 -5.75045075 0.92306020 1.4
2009 Q2 -0.45174364 0.70899368 -3.00372447 16.09059408 0.8
2009 Q3 0.60491332 -1.10343180 1.39880419 -24.49229966 0.3
2009 Q4 -0.01115014 -0.13213193 1.54400617 0.84829220 0.1
2010 Q1 0.53481740 0.10094986 1.88006931 -5.54399051 0.0
2010 Q2 0.81040406 1.29229259 2.05402479 11.65612884 -0.5
2010 Q3 0.64501881 0.49678098 1.42683671 -0.35208609 0.1
2010 Q4 1.01833874 0.69495229 0.37927209 -3.27335958 -0.2
2011 Q1 0.50041315 1.21571502 0.50174040 14.33860193 -0.3
2011 Q2 0.20141978 -0.15658108 0.21878696 -4.07705131 0.1
2011 Q3 0.43372599 0.52891255 1.01113866 2.72250400 -0.1
2011 Q4 0.33593895 0.06074719 0.85151692 -3.45447712 -0.5
2012 Q1 0.60108995 1.62204885 0.88651817 17.62530510 -0.3
2012 Q2 0.16942956 0.76689543 0.62923586 8.96949710 0.0
2012 Q3 0.26416034 -0.05071452 0.07880166 -3.04922177 -0.4
2012 Q4 0.27877186 2.59106697 0.63305509 29.04670355 0.1
2013 Q1 0.46861292 -4.26525047 0.67713243 -68.78826698 -0.4
2013 Q2 0.20545802 0.58146541 0.30744961 7.81647729 0.0
2013 Q3 0.46641787 0.58328912 0.23440888 3.49400682 -0.3
2013 Q4 0.83917367 0.21494896 0.79208722 -11.27661450 -0.5
2014 Q1 0.47345118 1.10369487 0.54709166 13.52020248 0.0
2014 Q2 0.93375698 1.29390492 1.33801074 8.24404770 -0.6
2014 Q3 0.91687178 0.99853396 0.62352731 2.46195256 -0.2
2014 Q4 1.12533250 1.04641801 0.90355427 -1.51305022 -0.3
2015 Q1 0.59624005 0.49040680 -0.46710878 -0.75840017 -0.2
2015 Q2 0.70814389 0.95495949 -0.69702162 5.02391773 -0.1
2015 Q3 0.66496956 0.80166267 0.38060610 3.18092976 -0.3
2015 Q4 0.56167978 0.74006260 -0.84554638 3.48278601 0.0
2016 Q1 0.40468216 0.51902540 -0.41793048 2.23653405 0.0
2016 Q2 1.04770741 0.72372078 -0.20331883 -2.72150106 -0.1
2016 Q3 0.72959779 0.64470081 0.47491844 -0.57285793 0.0
plot(uschange) # 'classical view'
autoplot(uschange) # 'all the series together'
# if we want to see just the first two series
par(mfrow=c(2,1))
plot(uschange[,1])
plot(uschange[,2])
tsdisplay(uschange[,1])
# to go back to 1 panel
par(mfrow=c(1,1))
##estimate an ARMAX model
armax1<- Arima(uschange[,1], xreg=uschange[,2], order=c(1,0,1))
res1<- residuals(armax1)
Acf(res1)
fitted(armax1)
Qtr1 Qtr2 Qtr3 Qtr4
1970 0.78006603 0.76998241 0.77997918 0.49470889
1971 0.74691458 0.95367578 0.76138513 0.87286648
1972 0.90802836 1.05689010 1.35176090 1.67290915
1973 1.16367516 1.23640518 0.78531105 0.79005897
1974 0.05003515 0.02887723 0.37542403 0.39653396
1975 0.01085034 1.17153692 0.24080247 0.93242178
1976 0.99640757 1.06214642 0.99472577 0.93136100
1977 0.85572449 1.09847602 0.97310973 1.02294306
1978 0.89498278 0.84279905 1.06793243 0.86730000
1979 0.87840579 0.42871159 0.56972244 0.68221667
1980 0.57340521 0.10300465 0.11068218 0.60940787
1981 0.43861630 0.52865939 0.79754365 0.41344059
1982 0.18783151 0.46954643 0.49131738 0.58362522
1983 0.95498317 0.89199407 1.24400153 1.39600580
1984 1.33190327 1.08208133 1.06911496 0.85908042
1985 0.74756500 1.31211568 0.74540820 1.23350168
1986 0.94718691 0.91606619 0.83474392 0.93429958
1987 0.92479494 0.36025112 1.11010775 0.99798195
1988 0.81430804 0.96157003 0.86476627 0.84641084
1989 0.94794685 0.53818285 0.70483443 0.79746028
1990 0.71661438 0.71355530 0.52369366 0.36479282
1991 0.31875176 0.32696302 0.46657727 0.60408691
1992 0.72976714 0.84119795 0.68756805 0.75856897
1993 0.87592653 0.66138154 0.68859617 0.87554555
1994 0.85504579 0.95994479 0.80156428 0.92224067
1995 0.86455278 0.58247944 0.78777694 0.75592063
1996 0.80440335 0.82728814 0.84656787 0.72326657
1997 0.82077487 0.83972711 0.79895963 1.04971381
1998 1.16710358 0.96916585 1.05446203 1.01986675
1999 1.07924160 0.86920134 1.03626172 1.13920607
2000 1.24577940 1.07983809 1.01776584 0.81608611
2001 0.89572256 0.59217313 0.96084875 0.19014878
2002 1.19322752 0.61764119 0.49594910 0.67128775
2003 0.61055054 0.79469592 0.91574479 0.81404964
2004 0.83743906 0.89386774 0.73441361 0.90861396
2005 0.51585232 0.87493303 0.84602049 0.85710258
2006 0.99228742 0.64174490 0.64139562 0.81848258
2007 0.75090681 0.61238757 0.57143827 0.52902716
2008 0.57538416 0.68538909 -0.17152326 0.32110519
2009 -0.08758067 0.15775900 -0.18758071 0.29165312
2010 0.32840865 0.65512112 0.58502461 0.65935684
2011 0.85147481 0.51911677 0.58684456 0.49275126
2012 0.77004410 0.60517415 0.38933485 0.88819759
2013 -0.45825324 0.69458103 0.59071556 0.52254602
2014 0.78259775 0.75487573 0.76063918 0.81156876
2015 0.77871811 0.80690024 0.75281487 0.72333875
2016 0.65205361 0.64574770 0.74192959
plot(uschange[,1])
lines(fitted(armax1), col=2)
armax2<- Arima(uschange[,1], xreg=uschange[,2], order=c(1,0,2))
res2<- residuals(armax2)
Acf(res2)
fitted(armax2)
Qtr1 Qtr2 Qtr3 Qtr4
1970 0.78312833 0.78018625 0.80479451 0.41930867
1971 0.84727727 0.77029485 0.84699950 0.95331366
1972 0.84462865 1.10359869 1.37654140 1.83317035
1973 1.13786583 1.20826089 0.98517520 0.70586675
1974 -0.05888372 -0.10652716 0.15252516 0.37252999
1975 0.26206496 1.00859262 0.13720674 0.91699737
1976 1.21948453 1.08521548 1.14267609 0.94879507
1977 0.77480045 1.08771984 1.08065070 0.98559686
1978 0.77381614 0.89892787 0.91640998 1.05941956
1979 0.89261857 0.38091069 0.53759910 0.54387460
1980 0.62276776 0.12597681 0.28963813 0.22350695
1981 0.29439649 0.66462244 1.02513172 0.44125915
1982 0.22393882 0.25312328 0.48707851 0.59336917
1983 0.90313806 1.08127964 1.25264967 1.50432058
1984 1.44036881 1.17436646 0.95189308 0.84432657
1985 0.62302160 1.25286288 0.88979506 1.13370425
1986 1.15901792 0.86946926 0.70947373 0.78744003
1987 1.06836399 0.45881327 0.90243930 1.03659542
1988 0.96270927 0.78370422 0.93047644 0.86924854
1989 0.91533693 0.58627680 0.66536135 0.72658394
1990 0.77608713 0.67685037 0.55014995 0.30925675
1991 0.39935061 0.27248164 0.26732620 0.62816783
1992 0.89742511 0.70072700 0.80056755 0.77549795
1993 0.88877001 0.77843395 0.62254092 0.81818587
1994 0.90212475 0.96475852 0.84927848 0.92193247
1995 0.82302821 0.63143690 0.68951971 0.73146966
1996 0.84367658 0.83076836 0.83169021 0.77671808
1997 0.80639086 0.80452612 0.86904546 0.93489503
1998 1.25190643 1.05908067 0.96434947 1.05049178
1999 1.09414605 0.93073127 0.96307717 1.17990917
2000 1.25691761 1.07448753 1.05545026 0.78640224
2001 0.84779060 0.60439534 0.96407640 0.13244226
2002 0.98024037 0.82330683 0.48567225 0.59802581
2003 0.62215476 0.83382918 0.88479540 0.77211433
2004 0.92453660 0.92649952 0.77051709 0.86179287
2005 0.46441540 0.86776267 0.86193275 0.90781283
2006 1.07577457 0.54136888 0.60248018 0.82077639
2007 0.72101375 0.64631554 0.58322422 0.49104710
2008 0.58814023 0.74866969 -0.28741575 0.27567649
2009 -0.01283404 0.01246222 -0.24998664 0.14688615
2010 0.47625237 0.73744143 0.63427916 0.71035984
2011 0.86308493 0.59824770 0.60961002 0.44522796
2012 0.76540887 0.60120651 0.39972191 0.89227477
2013 -0.52514212 0.53160833 0.74818031 0.56648441
2014 0.74603201 0.83062237 0.73814649 0.80864084
2015 0.76900278 0.88423189 0.77020729 0.70789375
2016 0.63785994 0.64367730 0.65987214
plot(uschange[,1])
lines(fitted(armax2), col=2)
AIC(armax2)
[1] 325.908
AIC(arima1)
[1] 353.1935
# Note: forecasts are available only if we have future values of 'income'
# procedure also available with auto.arima
auto.arima<- auto.arima(uschange[,1], xreg=uschange[,2])
Quarterly international arrivals (in thousands) to Australia from Japan, New Zealand, the UK and the US. 1981Q1 - 2012Q3
?arrivals
autoplot(arrivals)
autoplot(arrivals[,c(1,2)])
Japan<- arrivals[,1]
NZ<- arrivals[,2]
UK<- arrivals[,3]
US<- arrivals[,4]
# #we try with a simple arima model, ARMAX
auto.a<- auto.arima(NZ, xreg=Japan)
AIC(auto.a)
[1] 973.6105
# We try with a regression model with trend, season, and external variable 'Japan'
mod<- tslm(NZ~ trend+season+Japan)
summary(mod)
Call:
tslm(formula = NZ ~ trend + season + Japan)
Residuals:
Min 1Q Median 3Q Max
-40.32 -14.12 -2.51 12.70 55.63
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.73757 5.21347 0.909 0.36531
trend 2.21843 0.05539 40.052 < 2e-16 ***
season2 36.35773 5.01769 7.246 4.38e-11 ***
season3 63.33026 4.92212 12.866 < 2e-16 ***
season4 45.09142 4.95892 9.093 2.33e-15 ***
Japan -0.10039 0.03232 -3.106 0.00236 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.67 on 121 degrees of freedom
Multiple R-squared: 0.9482, Adjusted R-squared: 0.946
F-statistic: 442.9 on 5 and 121 DF, p-value: < 2.2e-16
fitted(mod)
Qtr1 Qtr2 Qtr3 Qtr4
1981 5.473922 44.596412 73.702533 56.744161
1982 14.111305 53.340009 82.418522 65.472097
1983 22.628218 62.051280 91.008119 74.004269
1984 31.251144 70.517295 99.486582 82.318292
1985 39.414981 78.991039 108.118945 90.566457
1986 47.510752 86.651414 116.267522 98.337764
1987 55.004376 94.934013 122.563578 104.736118
1988 61.010502 99.657338 127.933523 110.404225
1989 68.491577 108.845668 136.838850 120.604200
1990 74.807108 114.657135 141.688969 126.039901
1991 83.682418 122.719275 149.358178 132.036088
1992 87.828292 128.898776 156.225560 140.169807
1993 95.524004 136.766759 164.328660 147.883890
1994 103.102559 145.007697 171.680632 155.159866
1995 111.462962 152.265100 177.464491 163.072122
1996 117.239090 161.136997 185.808128 172.519362
1997 126.483740 170.719062 193.849186 181.070508
1998 137.181756 179.467177 205.411471 191.863193
1999 147.224618 189.712227 215.568177 201.293668
2000 156.278926 198.601893 225.039210 208.014379
2001 164.989795 207.314770 233.437962 222.448459
2002 175.508110 215.982270 242.616856 225.373172
2003 184.495357 230.450182 254.243793 234.589712
2004 192.461222 234.231033 260.581310 243.709174
2005 200.494450 244.901943 269.998834 253.592614
2006 209.805860 254.432307 280.136766 263.547634
2007 221.025108 265.056434 290.474578 274.694299
2008 232.838170 276.199485 302.070293 287.263910
2009 244.765880 287.395341 314.315139 297.607746
2010 253.343530 295.716391 320.226796 306.002282
2011 263.841867 306.382683 331.249779 315.888633
2012 272.184400 314.617597 339.578056
plot(NZ)
lines(fitted(mod), col=2)
plot(residuals(mod))
# analysis of residuals: autocorrelation?
tsdisplay(residuals(mod))
# fit an arima model to residuals
aar<- auto.arima(residuals(mod))
fitted(aar)
Qtr1 Qtr2 Qtr3 Qtr4
1981 43.6224119 42.8277172 12.1263286 5.1327012
1982 39.2568501 30.9438192 -7.0955060 -9.0387482
1983 19.5134000 9.3981650 -17.6940294 -17.4483315
1984 15.6074063 -2.5341788 -18.3538201 -19.9414017
1985 1.4237203 -4.5966746 -32.7036950 -33.3063718
1986 7.9409396 -6.3701574 -23.8781040 -8.2219509
1987 18.2048570 8.7177606 -10.6426164 6.1043757
1988 31.7536169 19.4098627 9.5556372 40.0743379
1989 42.7494458 21.4951003 5.9384598 -11.5850634
1990 9.7508942 -3.2840699 -16.1079337 -4.5609618
1991 -2.7173333 2.6010011 -10.8479635 -5.0272017
1992 8.0078896 -17.3846533 -18.6241736 -19.5752213
1993 -6.4796214 -9.0315699 -16.7577210 -4.0981064
1994 -5.3500015 -20.1121535 -24.2384935 -22.4297611
1995 -21.0249533 -27.2382296 -21.7057467 -12.6737140
1996 -5.7705269 -3.2515084 7.0836384 6.5702005
1997 13.2841639 1.9723677 0.9505590 -5.1161518
1998 11.7656134 -19.7140471 10.7822451 -2.7909927
1999 -1.5343071 -0.9442189 -1.6686489 -13.3979079
2000 -1.5534222 -12.9213959 15.9815234 -3.6481466
2001 10.8415074 12.1928739 4.5535895 -0.2197456
2002 -14.2595975 -13.9064412 -15.7578096 -27.9521114
2003 -11.7220191 -21.1521442 -22.1043363 -11.5919890
2004 -11.4794723 -5.6563040 15.6924561 37.1415765
2005 25.0658788 27.7020028 40.5272783 39.4642352
2006 13.4666275 8.5647731 33.0138668 15.9476795
2007 4.0567564 14.3061907 42.3182107 27.2034387
2008 0.0918427 15.4760632 17.3447321 -0.2983433
2009 -11.0766984 -0.1767021 -0.5943938 5.4166843
2010 -24.2731970 -5.7687911 -1.5124909 6.3192243
2011 -14.2285524 -7.6044283 4.0132082 -2.5252514
2012 -26.7724317 -9.4660669 -9.6347055
# complete the analysis by summing predictions made with linear model and arma on residuals
plot(NZ)
lines(fitted(mod)+fitted(aar), col=3)
# Notice the difference between the two methods: ARMAX and linear regression+ Arima on residuals
# #Another way of performing the same linear regression
tt<- (1:length(NZ))
seas <- factor(c(rep(1:4,length(NZ)/4),1:3)) #1:3 because there are three observations 'out'
mod2 <- lm(NZ~ tt+seas+Japan)
summary(mod2)
Call:
lm(formula = NZ ~ tt + seas + Japan)
Residuals:
Min 1Q Median 3Q Max
-40.32 -14.12 -2.51 12.70 55.63
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.73757 5.21347 0.909 0.36531
tt 2.21843 0.05539 40.052 < 2e-16 ***
seas2 36.35773 5.01769 7.246 4.38e-11 ***
seas3 63.33026 4.92212 12.866 < 2e-16 ***
seas4 45.09142 4.95892 9.093 2.33e-15 ***
Japan -0.10039 0.03232 -3.106 0.00236 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.67 on 121 degrees of freedom
Multiple R-squared: 0.9482, Adjusted R-squared: 0.946
F-statistic: 442.9 on 5 and 121 DF, p-value: < 2.2e-16
AIC(mod2)
[1] 1124.947
AIC(mod)
[1] 1124.947
Forecasting electricity demand
elecdaily
Time Series:
Start = c(1, 4)
End = c(53, 4)
Frequency = 7
Demand WorkDay Temperature
1.428571 174.8963 0 26.0
1.571429 188.5909 1 23.0
1.714286 188.9169 1 22.2
1.857143 173.8142 0 20.3
2.000000 169.5152 0 26.1
2.142857 195.7288 1 19.6
2.285714 199.9029 1 20.0
2.428571 205.3375 1 27.4
2.571429 228.0782 1 32.4
2.714286 258.5984 1 34.0
2.857143 201.7970 0 22.4
3.000000 187.6298 0 22.5
3.142857 254.6636 1 30.0
3.285714 322.2323 1 42.4
3.428571 343.9934 1 41.5
3.571429 347.6376 1 43.2
3.714286 332.9455 1 43.1
3.857143 219.7517 0 23.7
4.000000 186.9816 0 22.3
4.142857 228.4876 1 24.0
4.285714 222.3650 1 23.1
4.428571 220.4669 1 23.3
4.571429 240.5806 1 29.6
4.714286 224.6652 1 22.1
4.857143 181.1498 0 20.3
5.000000 182.2231 0 27.0
5.142857 230.2226 0 34.5
5.285714 309.3744 1 41.4
5.428571 243.2239 1 24.5
5.571429 258.3896 1 29.2
5.714286 259.3053 1 26.8
5.857143 241.3185 0 29.2
6.000000 264.5050 0 38.6
6.142857 264.9556 1 33.9
6.285714 221.6807 1 25.7
6.428571 237.1498 1 30.5
6.571429 276.4134 1 35.1
6.714286 283.3155 1 35.6
6.857143 272.5084 0 40.0
7.000000 243.6298 0 39.5
7.142857 223.0845 1 24.7
7.285714 238.7871 1 26.2
7.428571 252.2366 1 26.9
7.571429 250.9030 1 26.6
7.714286 242.4558 1 26.6
7.857143 214.2726 0 28.7
8.000000 175.2935 0 23.0
8.142857 216.4104 1 25.3
8.285714 232.6376 1 30.2
8.428571 222.3622 1 25.4
8.571429 210.2031 1 20.3
8.714286 211.6370 1 20.4
8.857143 186.7995 0 19.6
9.000000 178.2322 0 22.9
9.142857 221.6854 1 24.3
9.285714 237.5595 1 32.7
9.428571 220.5093 1 21.8
9.571429 216.7880 1 21.5
9.714286 214.8875 1 23.5
9.857143 193.5388 0 23.1
10.000000 183.3550 0 22.3
10.142857 228.3825 1 25.0
10.285714 255.3741 1 33.3
10.428571 233.9622 1 25.6
10.571429 216.4703 1 20.8
10.714286 215.9607 1 22.5
10.857143 194.5699 0 27.8
11.000000 200.7804 0 32.7
11.142857 213.9160 0 32.1
11.285714 240.5818 1 31.1
11.428571 220.3766 1 21.8
11.571429 219.1226 1 22.7
11.714286 223.9310 1 26.2
11.857143 191.9780 0 29.1
12.000000 172.9602 0 20.3
12.142857 212.6693 1 24.7
12.285714 216.7992 1 24.1
12.428571 215.9562 1 22.2
12.571429 226.8287 1 30.7
12.714286 221.9756 1 22.9
12.857143 186.2491 0 19.4
13.000000 179.2746 0 18.7
13.142857 213.5667 1 18.8
13.285714 214.6133 1 24.4
13.428571 217.4129 1 26.9
13.571429 221.7905 1 25.3
13.714286 214.8814 1 21.5
13.857143 185.2121 0 23.6
14.000000 181.7951 0 25.9
14.142857 230.6637 1 30.1
14.285714 253.2901 1 32.9
14.428571 235.5289 1 26.0
14.571429 222.7004 1 18.2
14.714286 218.9320 1 19.7
14.857143 192.1518 0 24.1
15.000000 182.5951 0 24.3
15.142857 218.4745 1 24.5
15.285714 226.2018 1 22.0
15.428571 225.5993 1 17.8
15.571429 223.8808 1 18.9
15.714286 220.3987 1 18.9
15.857143 189.7302 0 22.4
16.000000 180.7974 0 21.5
16.142857 212.9014 1 20.7
16.285714 216.1034 1 21.1
16.428571 218.3603 1 21.6
16.571429 215.1137 1 24.3
16.714286 177.7474 0 18.8
16.857143 183.7940 0 17.3
17.000000 178.9498 0 18.3
17.142857 181.3320 0 18.1
17.285714 219.8456 1 17.6
17.428571 218.9037 1 19.2
17.571429 217.9708 1 18.3
17.714286 189.6816 0 20.5
17.857143 189.3522 0 19.2
18.000000 186.1761 0 17.2
18.142857 222.7253 1 21.3
18.285714 226.8940 1 22.3
18.428571 227.8383 1 17.6
18.571429 234.1514 1 15.6
18.714286 234.8755 1 13.9
18.857143 205.5938 0 17.0
19.000000 203.6670 0 13.7
19.142857 237.8857 1 16.4
19.285714 240.6153 1 16.0
19.428571 240.3634 1 16.2
19.571429 241.3502 1 18.4
19.714286 236.7497 1 18.3
19.857143 206.3455 0 16.3
20.000000 192.7160 0 19.7
20.142857 229.3544 1 19.9
20.285714 230.0164 1 20.1
20.428571 224.4637 1 20.2
20.571429 216.9592 1 22.4
20.714286 212.6305 1 22.8
20.857143 188.9991 0 20.0
21.000000 181.5075 0 21.4
21.142857 217.0623 1 22.8
21.285714 224.3771 1 21.6
21.428571 221.9410 1 23.1
21.571429 222.4899 1 20.9
21.714286 221.5659 1 20.3
21.857143 195.6669 0 17.1
22.000000 190.3111 0 20.2
22.142857 226.0510 1 20.8
22.285714 226.9175 1 21.7
22.428571 230.2927 1 18.0
22.571429 234.6604 1 15.8
22.714286 231.8486 1 17.7
22.857143 201.0383 0 18.4
23.000000 199.7129 0 14.9
23.142857 232.5876 1 16.4
23.285714 235.0114 1 16.4
23.428571 230.9970 1 18.1
23.571429 232.0875 1 16.7
23.714286 234.6715 1 15.4
23.857143 202.4323 0 16.8
24.000000 196.7197 0 14.9
24.142857 206.5046 0 15.4
24.285714 235.7782 1 18.3
24.428571 241.0994 1 16.2
24.571429 243.6167 1 16.2
24.714286 237.2126 1 16.5
24.857143 212.1366 0 15.2
25.000000 203.3582 0 16.0
25.142857 241.3206 1 17.0
25.285714 244.7568 1 15.8
25.428571 244.5883 1 14.2
25.571429 252.8974 1 12.2
25.714286 239.6948 1 15.1
25.857143 211.0494 0 15.0
26.000000 205.6888 0 15.4
26.142857 245.3576 1 14.4
26.285714 261.7745 1 14.3
26.428571 248.2476 1 15.7
26.571429 243.3896 1 15.6
26.714286 241.5915 1 14.5
26.857143 218.8741 0 12.7
27.000000 220.2953 0 11.0
27.142857 255.0056 1 13.1
27.285714 254.8101 1 13.1
27.428571 245.0865 1 16.6
27.571429 247.8669 1 12.4
27.714286 249.0520 1 12.4
27.857143 218.4067 0 15.2
28.000000 212.9148 0 13.9
28.142857 241.2600 1 15.7
28.285714 242.9725 1 14.3
28.428571 257.5991 1 13.3
28.571429 257.9806 1 14.7
28.714286 247.6547 1 14.4
28.857143 227.0290 0 11.1
29.000000 220.9397 0 14.9
29.142857 258.8287 1 13.1
29.285714 264.1827 1 12.9
29.428571 253.5232 1 14.6
29.571429 259.8050 1 12.1
29.714286 256.2355 1 13.8
29.857143 223.4825 0 13.2
30.000000 215.7748 0 13.2
30.142857 258.8278 1 12.9
30.285714 268.9846 1 12.2
30.428571 267.2355 1 13.4
30.571429 269.3544 1 14.9
30.714286 253.2748 1 15.9
30.857143 223.1568 0 14.9
31.000000 211.2160 0 16.8
31.142857 242.0725 1 16.7
31.285714 246.4046 1 16.5
31.428571 237.8841 1 18.9
31.571429 239.6187 1 19.3
31.714286 263.3148 1 9.8
31.857143 227.5769 0 13.4
32.000000 221.6350 0 14.2
32.142857 256.6197 1 15.2
32.285714 250.4816 1 16.5
32.428571 246.8155 1 15.9
32.571429 252.0098 1 13.6
32.714286 243.5690 1 15.9
32.857143 217.1874 0 14.6
33.000000 212.9136 0 13.6
33.142857 261.0586 1 12.7
33.285714 262.4751 1 12.7
33.428571 253.5156 1 14.5
33.571429 249.8530 1 14.8
33.714286 244.9549 1 16.9
33.857143 221.1463 0 13.3
34.000000 210.5977 0 15.7
34.142857 250.2505 1 13.6
34.285714 248.2650 1 13.3
34.428571 247.3034 1 13.3
34.571429 242.7937 1 17.4
34.714286 231.1652 1 20.6
34.857143 202.5800 0 19.9
35.000000 200.8607 0 16.3
35.142857 239.1268 1 16.1
35.285714 236.9323 1 18.7
35.428571 233.0740 1 19.5
35.571429 233.1221 1 20.0
35.714286 233.8208 1 15.0
35.857143 193.0105 0 21.0
36.000000 189.3291 0 21.0
36.142857 236.2696 1 17.8
36.285714 238.2240 1 16.1
36.428571 241.4517 1 14.8
36.571429 241.7054 1 15.9
36.714286 233.5166 1 16.7
36.857143 198.1050 0 18.9
37.000000 186.0104 0 20.7
37.142857 218.2276 1 20.8
37.285714 222.7274 1 20.7
37.428571 222.0213 1 19.1
37.571429 226.7373 1 17.1
37.714286 231.3683 1 15.9
37.857143 196.5834 0 18.5
38.000000 185.0438 0 20.1
38.142857 225.4641 1 18.0
38.285714 235.6982 1 14.8
38.428571 236.1329 1 16.2
38.571429 241.1509 1 13.7
38.714286 230.3640 1 15.5
38.857143 199.8989 0 16.8
39.000000 187.0570 0 17.2
39.142857 217.4684 1 23.0
39.285714 215.1564 1 25.3
39.428571 215.4661 1 23.5
39.571429 217.7978 1 17.0
39.714286 214.7996 1 19.3
39.857143 184.7269 0 23.3
40.000000 177.3159 0 27.9
40.142857 207.7305 1 22.9
40.285714 218.1746 1 23.8
40.428571 224.6942 1 15.9
40.571429 224.9024 1 17.4
40.714286 212.3656 1 18.1
40.857143 183.7487 0 26.2
41.000000 173.5839 0 19.2
41.142857 221.1614 1 25.0
41.285714 216.7209 1 17.7
41.428571 221.2251 1 18.7
41.571429 217.7267 1 23.1
41.714286 213.9397 1 24.9
41.857143 186.2191 0 21.2
42.000000 181.4456 0 27.5
42.142857 222.1222 1 15.5
42.285714 228.7036 1 16.9
42.428571 225.6833 1 17.1
42.571429 224.5799 1 16.7
42.714286 214.5551 1 18.1
42.857143 191.7391 0 24.2
43.000000 183.9389 0 28.8
43.142857 213.2147 1 20.6
43.285714 222.9096 1 28.2
43.428571 233.0016 1 31.1
43.571429 224.6348 1 22.7
43.714286 227.9454 1 30.4
43.857143 193.2198 0 19.8
44.000000 182.9607 0 24.5
44.142857 212.9383 1 22.5
44.285714 220.9596 1 19.7
44.428571 218.7275 1 22.5
44.571429 218.7459 1 21.8
44.714286 226.7615 1 32.7
44.857143 187.8364 0 20.3
45.000000 186.3019 0 16.6
45.142857 200.5634 1 23.3
45.285714 186.9384 0 28.9
45.428571 206.6377 1 17.8
45.571429 211.9857 1 23.5
45.714286 227.2132 1 32.8
45.857143 205.4283 0 34.6
46.000000 176.0597 0 19.4
46.142857 211.9239 1 20.2
46.285714 215.3201 1 18.8
46.428571 222.2635 1 24.3
46.571429 239.6066 1 30.6
46.714286 223.6820 1 27.7
46.857143 189.5637 0 18.0
47.000000 182.0131 0 18.1
47.142857 208.3546 1 20.1
47.285714 211.8105 1 20.9
47.428571 217.4702 1 29.6
47.571429 228.7774 1 28.1
47.714286 215.2890 1 23.0
47.857143 196.0592 0 30.8
48.000000 187.2593 0 25.0
48.142857 222.5512 1 19.6
48.285714 213.9397 1 21.8
48.428571 218.0656 1 18.7
48.571429 214.9400 1 19.4
48.714286 212.5017 1 20.5
48.857143 193.7848 0 24.3
[ reached getOption("max.print") -- omitted 32 rows ]
plot(elecdaily)
# fit a quadratic regression model with ARMA errors using the auto.arima function
xreg <- cbind(MaxTemp = elecdaily[, 'Temperature'],
MaxTempSq = elecdaily[, 'Temperature']^2,
Workday = elecdaily[, 'WorkDay'])
fit <- auto.arima(elecdaily[, 'Demand'], xreg = xreg)
checkresiduals(fit)
Ljung-Box test
data: Residuals from Regression with ARIMA(2,1,2)(2,0,0)[7] errors
Q* = 28.229, df = 8, p-value = 0.0004326
Model df: 6. Total lags used: 14
# The model has some significant autocorrelation in the residuals
# using the estimated model we forecast 14 days ahead, setting for the next 14 days a temperature at a constant level of 26.
fcast <- forecast(fit,
xreg = cbind(MaxTemp=rep(26,14), MaxTempSq=rep(26^2,14),
Workday=c(0,1,0,0,1,1,1,1,1,0,0,1,1,1)))
plot(fcast)
Lagged Predictors
plot(insurance, main='Insurance advertising and quotations', xlab='year', xaxt='n')
axis(1, at=c(2002,2003,2004,2005),labels=c('2002','2003','2004','2005'), line=-1, lwd=0, lwd.ticks=0.5)
# we consider 1 lagged predictor
advert<- cbind(insurance[,2], c(NA,insurance[1:39,2]))
colnames(advert)<- paste('AdLag',0:1, sep='')
fit<- auto.arima(insurance[,1], xreg=advert[,1:2],d=0)
summary(fit)
Series: insurance[, 1]
Regression with ARIMA(3,0,0) errors
Coefficients:
ar1 ar2 ar3 intercept AdLag0 AdLag1
1.4117 -0.9317 0.3591 2.0393 1.2564 0.1625
s.e. 0.1698 0.2545 0.1592 0.9931 0.0667 0.0591
sigma^2 = 0.2165: log likelihood = -23.89
AIC=61.78 AICc=65.4 BIC=73.43
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.0126932 0.4344493 0.3501642 -0.1590732 2.611605 0.1020281 0.001076967
# we provide the forecast by assuming that advertising is 8 units in each future month
f<-forecast(fit,h=20, xreg=cbind(AdLag0=rep(8,20),AdLag1=c(advert[40,1],rep(8,19))))
plot(f)
Twitter Revenues
twitter<- read_excel('twitter.xlsx')
length(twitter$twitter)
[1] 46
tw<- (twitter$twitter)
# GGM
GGM_tw<- GGM(tw, prelimestimates=c(4.463368e+04, 0.001, 0.01, 1.923560e-03, 9.142022e-02))
Warning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs produced
summary(GGM_tw)
Call: ( Guseo Guidolin Model )
GGM(series = tw, prelimestimates = c(44633.68, 0.001, 0.01, 0.00192356,
0.09142022))
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-257.21 -64.74 20.95 9.17 70.88 265.11
Coefficients:
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error 103.2668 on 41 degrees of freedom
Multiple R-squared: 0.999862 Residual sum of squares: 437225.1
pred_GGM_tw<- predict(GGM_tw, newx=c(1:60))
pred_GGM_tw.inst<- make.instantaneous(pred_GGM_tw)
plot(tw, type= 'b',xlab='Quarter', ylab='Quarterly revenues', pch=16, lty=3, cex=0.6, xlim=c(1,60))
lines(pred_GGM_tw.inst, lwd=2, col=2)
lines(pred_GGM_tw.inst, lwd=2, col=3)
# Analysis of residuals
res_GGMtw<- residuals(GGM_tw)
acf<- acf(residuals(GGM_tw))
fit_GGMtw<- fitted(GGM_tw) ##they are cumulative values
fit_GGMtw_inst<- make.instantaneous(fit_GGMtw)
SARMAX refinement
# SARMAX model with external covariate 'fit_GGM'
s2 <- Arima(cumsum(tw), order = c(3,0,1), seasonal=list(order=c(3,0,1), period=4),xreg = fit_GGMtw)
summary(s2)
Series: cumsum(tw)
Regression with ARIMA(3,0,1)(3,0,1)[4] errors
Coefficients:
ar1 ar2 ar3 ma1 sar1 sar2 sar3 sma1 intercept xreg
1.8294 -1.2030 0.1683 -1.000 0.5144 0.6488 -0.2487 0.0400 21.0796 0.9983
s.e. 0.1519 0.2605 0.1508 0.027 0.5208 0.2650 0.3104 0.5359 10.4055 0.0010
sigma^2 = 2329: log likelihood = -242.57
AIC=507.14 AICc=514.91 BIC=527.26
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.2324242 42.69616 30.24614 0.6453786 2.042453 0.04979198 0.03877914
pres2 <- make.instantaneous(fitted(s2))
plot(twitter$twitter, type= 'b',xlab='Quarter', ylab='Quarterly revenues', pch=16, lty=3, xaxt='n', cex=0.6)
axis(1, at=c(1,10,19,28,37,46), labels=twitter$quarter[c(1,10,19,28,37,46)])
lines(fit_GGMtw_inst, lwd=1, lty=2)
lines(pres2, lty=1,lwd=1, col=2)
Interest in Zoom according to searches in Google ####
interest<- read.csv('zoomgoogle.csv')
int<- interest$zoom[90:180]
week <- as.Date(interest$week[90:180])
# we set a 'timing' more refined for graphical purposes only
tfine150 <- seq(1,150,0.01)
# Exploratory plots
plot(int, type= 'b',xlab='Week', ylab='Weekly Google searches', pch=16, lty=3, cex=0.6, xaxt='n')
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))
plot(cumsum(int), type= 'b',xlab='Week', ylab='Cumulative Google searches', pch=16, lty=3, cex=0.6, xaxt='n')
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))
# GBM with one exponential shock
gbme1.go<- GBM(int,shock = 'exp',nshock = 1,prelimestimates = c(4.046997e+03, 0.001, 0.1, 30,-01,5))
# Prediction with GBMe1
pred.gbme1go<- predict(gbme1.go, newx=tfine150)
pred.gbme1goi<- make.instantaneous(pred.gbme1go)[-1]
# Plots observed vs predicted
plot(int, type= 'b',xlab='Week', ylab='Weekly Google searches', pch=16, lty=3, cex=0.6, xaxt='n', col=colors[3])
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))
lines(tfine150[-1], pred.gbme1goi*100, lwd=2, col=1)
plot(cumsum(int), type= 'b',xlab='Week', ylab='Cumulative Google searches', pch=16, lty=3, cex=0.6, xaxt='n',col=colors[3])
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))
lines(tfine150, pred.gbme1go, lwd=2, col=1)
SARMAX refinement
fit.google<- fitted(gbme1.go)
s2 <- Arima(cumsum(int), order = c(1,0,1), seasonal=list(order=c(0,0,1), period=52), xreg = fit.google)
summary(s2)
Series: cumsum(int)
Regression with ARIMA(1,0,1)(0,0,1)[52] errors
Coefficients:
ar1 ma1 sma1 intercept xreg
0.9155 0.5564 -0.1156 -0.1915 1.0019
s.e. 0.0402 0.0782 0.2036 10.9422 0.0071
sigma^2 = 16.99: log likelihood = -257.29
AIC=526.58 AICc=527.58 BIC=541.65
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.02828682 4.007078 2.47122 0.1528441 1.082046 0.08839818 0.1306401
pres2 <- make.instantaneous(fitted(s2))
# Plots observed vs predicted with SARMAX refinement
plot(int, type= 'b',xlab='week', ylab='Weekly Google searches', pch=16, lty=3, cex=0.6, xaxt='n', col=colors[3])
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))
lines(tfine150[-1], pred.gbme1goi*100, lwd=1, lty=2)
lines(pres2, lty=1,lwd=1)
legend('topright', legend=c('GBMe1','GBMe1+SARMAX'), lty=c(2,1))
plot(cumsum(int), type= 'b',xlab='week', ylab='Cumulative Google searches', pch=16, lty=3, cex=0.6, xaxt='n', col=colors[3])
axis(1, at=c(1,18,36,54,72,90), labels=format(week[c(1,18,36,54,72,90)], '%d/%m/%y'))
lines(tfine150, pred.gbme1go, lwd=1, lty=2)
lines(cumsum(pres2), lty=1,lwd=1)
legend('bottomright', legend=c('GBMe1','GBMe1+SARMAX'), lty=c(2,1))